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Abstract. We present a time-splitting spectral scheme for the Maxwcll-Dirac system and similar time- 
splitting methods for the corresponding asymptotic problems in the semi-classical and the non-relativistic 
regimes. The scheme for the Maxwell-Dirac system conserves the Lorentz gauge condition, is unconditionally 
stable and highly efficient as our numerical examples show. In particular we focus in our examples on the 
creation of positronic modes in the semi-classical regime and on the electron-positron interaction in the 
non-relativistic regime. Furthermore, in the non-relativistic regime, our numerical method exhibits uniform 
convergence in the small parameter 8, which is the ratio of the characteristic speed and the speed of light. 
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1. Introduction and asymptotic scaling 

The Maxwell-Dirac system (MD) describes the time-evolution of fast, i.e. relativistic spin-1/2 particles, say 
electrons and positrons, within external and self-consistent electromagnetic fields. In Lorentz gauge it is 
given by the following set of equations: 

3 



(1.1) 



mi> = E ak (t 9 * ~ q{Ak + AT) ) ^+^ v + yex ^ + mc2 w> 



subject to Cauchy initial data: 

(V\ t=Q = V^(x), ^V-| t=0 = V (1) ( x )' 

lA| t=0 = A<°>(x), d t A| t=f) = A«(x), ^| t=0 ^ (0) (x). 

The particle- and current- densities p and J = (.71 , , J3) are defined by: 

(1.3) p:=q\ip\ 2 , j k :=qc(iP,a k ip) c *=qcip-a k xP, k = 1,2, 3, 
where the spinor field tp = ip(t,x) = (ipi, ipi, tps, ip4) T <G C 4 is normalized s.t. 

(1.4) / m, X )\ 2 d X = i, 
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with t, x = (xi,X2,X3), denoting the time - resp. spatial coordinates. Further, V(t, x) and V ex (x) € K 
are the self- consistent resp. external electric potential and j4fc(i,x) £ R, resp. ^^(x) € K, represents the 
/cth-components of the self-consistent, resp. external, magnetic potential, i.e. A = (Ai, A2, A3). Here and in 
the following we shall only consider static external fields. The complex- valued, Hermitian Dirac matrices, 
i.e. f3,a k , are explicitly given by: 

"-ft -l) ■ ■" 

with I2, the 2x2 identity matrix and o~ k the 2x2 Pauli matrices, i.e 

m> - 1 : ; 1 . °* ■- (° ;) ■ °> - 






Finally, the physical constants, appearing in ( 1.1 )-( 1.3 ) , are the normalized Planck's constant H, the speed 
of light c, the permittivity of the vacuum e , the particle mass m and its charge q. 

Additionally to ( |1.1| , we impose the Lorentz gauge condition 

(1.7) 3tV(t,x) + cdivA(i,x) = 0, 

for the initial potentials (x), (x), and A (0) (x), A (1) (x). That means 

-VW(x) + V ■ A(°>(x) = 0, AyC'fx) + -^-|V> (0) | 2 + -V ■ A«(x) = 0. 

C 47T60 c 

Then the gauge is henceforth conserved during the time-evolution. This ensures that the corresponding 
electromagnetic fields E, B are uniquely determined by 

(1.8) E(i,x) := -~dtA(t,x) - VV(*,x), B(i,x) := curl A(i,x). 

c 

Also it is easily seen that multiplying the Dirac equation with ip implies the following conservation law 

(1.9) dtp + div J = 0. 

The MD equations are the underlying field equations of relativistic quantum electro-dynamics, cf. |22) . where 
one considers the system within the formalism of second quantization. Nevertheless, in order to obtain a 
deeper understanding for the interaction of matter and radiation, there is a growing interest in the MD 
system also for classical fields, since one can expect at least qualitative results, cf. [T2]. Analytical results 



concerning local and global well-posedness of (1.1|-(1.3|, have been obtained in [TOl [TTJ [15l [16] . Also the 



rigorous study of asymptotic descriptions for the MD system has been a field of recent research. In particular 
the non- relativistic limit and the semi-classical asymptotic behavior (in the weakly coupled regime) have 
been discussed in 8, 215]. For the former case a numerical study can be found in [3]. Since our numerical 
simulations shall deal with both asymptotic regimes, let us discuss now more precisely the corresponding 
scaling for these physical situations. 

1.1. The MD system in the (weakly coupled) semi-classical regime. First, we consider the semi- 
classical or high-frequency regime of fast (relativistic) particles, i.e. particles which have a reference speed 
u«c. (Of course for particles with mass to > we always have v < c.) To do so we rewrite the MD system 
in dimensionless form, such that there remains only one positive real parameter 

(1.10) Kq = — . 
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As described in [23] , we obtain the following rescaled MD system: 

iK d t ^ = -ikq a ■ V^> - a • (A + A ex )ip + {V + V ex )tp + /3ip, 
il.ll] {(d tt -A)V = p, 

(d tt ~ A) A = J, 

where from now on we shall also use the shorthand notation a. ■ V := dk- Notice that if |g| = e, i.e. 
in the case of electrons or positrons where q equals the elementary charge ±e, the parameter kq 137 is 
nothing but the reciprocal of the famous fine structure constant. Thus for fast (relativistic) particles which 
are not too heavily charged, Kq in general is not small and therefore asymptotic expansions as k — > do 
not make sense. In order to describe the semi-classical regime we therefore suppose that the given external 
electromagnetic potentials are slowly varying w.r.t. the microscopic scales, i.e. V ext — V ext (x.e/Ko) and 
likewise A ext = A ext (xe / ko) , where from now on < e « 1 denotes the small semi-classical parameter. 
Here we fix kq and include it in the scaling which conveniently eliminates this factor from the resulting 
equations. Finally, observing the time-evolution on macroscopic scales we are led to 

(1.12) x=— x, t=—t. 

Kq Kq 

and we set 

(1.13) <A £ (i,xH(£) %(t^,x^) = (£) 3 %(*,x), 



in order to satisfy the normalization condition (1.4). Plugging this into |lTT]) and omitting all "~" we obtain 



the following semi-classically scaled MD system: 

!ie d t ip e = -iea-Vip E a- (A e + A ex ) ^ + (V E + V ex )ip e + pip e , 
(d tt ~A)V £ =e\r\ 2 , 
(d tt - A)A| = e(r, a k r)c^ k = 1, 2, 3, 

with < £ « 1. Note the additional factor e in the source terms appearing on the right hand side of the wave 
equations governing V 6 and A £ , which implies that we are dealing with a weak nonlinearity in the sense of 
[T31 121) . The scaled particle-density in this case is p £ := \ip £ \ 2 and we also have J e := ((tp e , a fe, e )c 4 )fe=i 2,3- 

Remark 1.1. Note that, equivalently, we could consider small asymptotic solutions ip e ~ 0(y/s) which 



again satisfy the semi-classical scaled MD system (1.14) but with source terms of order 0(1) in the wave 
equations. This point of view is adopted in |23j . 

1.2. The MD system in the non-relativistic regime. We shall also deal with the non-relativistic regime 
for the MD system, i.e. we consider particles which have a reference speed v <C c. Introducing a reference 
length L, time T and writing v — L/T, we rescale the time and the spatial coordinates in ( |1.1[ ) by 

(1-15) i=j, t = |. 



Moreover we set ip(t, x) = L 3 / 2 ip(t,Tt), such that (1.4) is satisfied, and we also rescale the electromagnetic 
potentials by 

(1.16) A (CT) (f,x) = AA (ea;) (i,x), y( ei '(a) = AF (ei) (i,x), 
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where A = q/ (AttLeq), cf. [21 16]. In this case we have again two important dimensionless parameters, namely 

(1.17) 5=-<l, k= 5-^. 

c qr 

Note that for w sa c we get k~hq. Choosing for convenience v = q 2 / (Airheo) and L = q/Aireo, we shall from 
now on denote by ip s (t, x) the rescaled wave function t/>(£, x), which is obtained for this particular choice of 
f = L/T. Then, similarly as before, ip s satisfies a dimensionless one-parameter model (again omitting all 
"~"), given by 



(1.18) 



id t ^ = -i a • W 4 - a • (A 5 + A ea; )V/ + (V 5 + V^)V 4 + ^/3^ 5 , 
(6 2 d tt -A)V s = \i> s \ 2 , 
[ (S 2 d tt - A)A 5 k = (V 5 , qV)c> fc = 1, 2, 3. 



In analogy to the semi-classical case, this system will be called the non-relativistically scaled MD system. 
In this case the scaled particle density is p 5 := \ip 5 \ 2 , whereas 3 s := <y _1 ((^>*, a fe- <5 )c 4 )fc=i.2,3. Note that in 
this scaling J 5 ~ O(l), A 5 ~ C((5), (due to a rather complex cancellation mechanism already known in the 
linear case cf. [5]) such that the magnetic field is a relativistic effect which does not appear in the zeroth 
order approximation of the MD system, cf. [SJ |H1 EI] (see also [7] for a similar study) . 

As in the corresponding numerical simulations for semi-classical nonlinear Schrodinger equations, cf. [T] , the 
main difficulty is to find an efficient and convenient numerical scheme with best possible properties in the 
limiting regimes e — > and 8 — > 0, i.e. in particular with uniform convergence properties in 8. 

In the following we present a time-splitting spectral method for the MD system, and its semi-classical and 
non-relativistic limiting systems. The time-splitting spectral methods have been proved to be the best 
numerical approach to solve linear and nonlinear Schrodinger type systems in the semi-classical regime, cf. 
[TJ[2]. Besides the usual properties of the time-splitting spectral method, such as the conservation of the 
Lorentz gauge condition and the unconditionally stability property, here we shall pay special attention to 
its performance in both the semi-classical and non-relativistic regimes. Note that in particular the semi- 
classical asymptotics has not been studied in [3]. The method proposed here is similar to the one used for 
the Zakharov system in [T5] . A distinguished feature of the scheme developed in [TH] is that it can be used, in 
the sub-sonic regime, with mesh size and time step independent of the subsonic parameter, a possibility not 
shared by works before [H [S] . For the MD system, our time splitting spectral method allows the use of mesh 
size and time steps independent of the relativistic parameter 8, allowing coarse grid computations in this 
asymptotic regime. This is achieved by the Crank-Nicolson time discretization for the Maxwell equations, a 
scheme shown to perform better for wave equations in the subsonic regime than the exact time integration, 
as studied in [TB]. For the same reason, the previously proposed time-splitting spectral method for the MD 
system in [3] does not possess this property since it uses the exact time integration for the Maxwell equations. 

The paper is now organized as follows: In section [2] we give the time-splitting spectral method for the MD 
system and one simple example to show the reliability, efficiency and the convergent rate of our method. Our 
method has spectral convergence for space discretization and second order convergence for time discretization. 
In section[3]and|4j we discuss the time-splitting methods for the asymptotic systems (the semi-classical regime 
and non-relativistic regime) and give some examples for them respectively. We conclude the paper in section 
5. 
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2. A TIME-SPLITTING SPECTRAL METHOD FOR THE MAXWELL-DlRAC SYSTEM 

2.1. A time-splitting method. Before we describe our time-splitting spectral method, we combine the 
rescaled MD system (1.14 1 and (1.18), using two parameters: 

f IF I 

is d t i> = a ■ W - a • (A e + A ex )^ + (V £ + V ex )^ + -W, 
( 2A ) \ (5 2 d ft -A)V = e\iP\ 2 , 

k (S 2 d tt - A)A k = e(iP, q^) C 4, k = 1, 2, 3. 

In the following we shall denote by 

(2.2) 2?aW : = » ' H V - Ae "( x ))^ + W + V ex {yi)i), 

the standard Dirac operator with (external) electromagnetic fields, D := — iV. The corresponding 4x4 
matrix- valued symbol is given by 

(2.3) X> A (0 = a • (C - A e *(x)) + ^ + V^x), 
where x, £ eK 3 . Likewise the free Dirac operator will be written as 

(2.4) V (D x )^:=-ia-V^ + ^ 
Its symbol admits a simple orthogonal decomposition given by 

(2-5) Z>„(0 = a • £ + p = A (On+(0 - VOHq (0, 

where 

(2-6) A (0 := Vllf+l, 

and 



(2.7) 



1 



Lt ± 



1 



Ao(0 



Z>o(0 



The time-splitting scheme we propose is then as follows 
Step 1. Solve the system 

1 



ieft^ - ^X» (fei? x )^ = 0, 
(5 2 d tt -A)V = s\^\ 2 , 
y (S 2 d tt - A)A k = e (V, «V) , * = 1, 2, 3, 



on a fixed time-interval At, using the spectral decomposition (2.5). 
Step 2. Then, in a second step we solve 



(2.9) 



ied t ip + a ■ (A + A ex )V - (V + V ex )4> = 0, 



on the same time-interval, where the solution obtained in step 1 serves as initial condition for step 2. Also 
the fields A, V are taken from step 1. It is then easy to see that this scheme conserves the particle density 
and the Lorentz gauge. 
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2.2. The numerical algorithm. In the following, for the convenience of computation, we shall deal with 
the system (2.1 1 on a bounded domain, for example, on the cubic domain 

(2.10) n = { X =(x u x 2 ,x 3 ) | aj <x 3 <bj,j = 1,2,3}, 

imposing periodic boundary conditions. We choose the time step At = T/M and spatial mesh size Axj — 
(bj — aj)/Nj, j = 1,2,3, in a^-direction, with given M,Nj € N and [0,T] denoting the computational time 
interval. Further we denote the time grid points by 

1 ' 



(2.11) 



t n = nAt, t n+ i/ 2 = ( n 



At, t = 0,l,. 



and the spatial grid points by 

(2.12) x m = (xi, TOl , x 2 ,mz, x 3>m3 ), where Xj, mj ■= a? +mjAxj, j = 1,2,3, 
and m = (7711,7712,7713) £ M., with 

(2.13) M = {(Tni,m 2 ,m 3 ) < mj < Nj, j = 1, 2, 3 j . 

In the following let Vl/^, V£> and AJ^ be the numerical approximations of ip(t n ,yi m ), V(t n ,x m ), and 
A(t„,x m ), respectively. Suppose that we are given $ n , V n , and A", then we obtain 9 n+1 , V n+1 and 
A n+1 as follows: 

Step 1. For the first step we denote the value of ^ at time t by $(£). Then we approximate the spatial 
derivative in (|2.8|) by the spectral differential operator. More precisely we first take a discrete Fourier 

9 t l> = — % - (e5a ■ £ + p) $ = Mi$, 

(6 2 d tt + \tf)V = efip, 

_ (<5 2 d tt + |£| 2 )i fe =eM§), forfc = l,2,3, 

where / is the DFT of function /. As the matrix Mi 6 C 4x4 j s diagonalizable, i.e. there exists a Hermitian 
matrix Z?i such that 

(2.15) D^Mi-Di = diag [A, A, —A, —A] = A 

is a purely imaginary diagonal matrix with entries 

(2.16) 

Then the value of <E> at time t n+ \ is given by 
= D x exp (AAi) L»f * n 



transform (DFT) of (2.8 1: 



(2.14) 



(2.17) 

where 
(2.18) 



/ c A - is\ -«^sa6 

ca-7s a sSs\(( 2 -i£i) 

-ie5s\£ 3 -£<5s(6 + «£i) c A + 7s A 



c A := cos(-iAAt), s A := sin(-iAA£)(l + \e6£\*) 



7£(5sa6 



ca + is A / 

2\-l/2 
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Then we obtain the value of $ n+1 by an inverse discrete Fourier transform (IDFT). Hence from (2.14), we 
can find the values of V and A by the Crank-Nicolson scheme, i.e. 



(2.19) 

and 

(2.20) 





1 



1 - 



At|g| 2 



At 2 |C| 2 

4(5 2 
At|g| 2 
S 2 




1 




n+l 



where for k = 1, 2, 3, we denote 

(2.21) p n = \^ n \ 2 , p n+1 = |$ n+1 | 2 , J£ = 5- 1 (y n ,a k V n ), 3 n k +1 = S~ 1 ($ n+1 ,a k $ n+1 ) . 
Performing an IDFT of V n+1 and A"+\ we finally obtain V n+1 and A n+1 . 

Step 2. Since V and do not change in Step 2, we only have to update "J. First we shall rewrite the 
equation (2.9| in the following form: 

(2.22) 9 t * = - a ■ (A + A 6 *)* - (V + V^)* = M 2 #. 
Then there exists again a Hermitian matrix £) 2 such that 

(2.23) D%M 2 D 2 = diag Ml , /i 2 , /x 2 ] = 6, 
where is a purely imaginary diagonal matrix with 

(2.24) ^^~ l -((V + V ex )-\A + A ex \), m = 
Hence, the value of \1/ at time t n+ \ is given by 



((V + V ex ) + I A + A ex |) 



n+l 



D 2 exp (QAt) D 2 r & n+1 

/ C1+C2— i(«l+fl2) 



(2.25) 







ci+c 2 — i(^i+s 2 ) 
2 

Ai-iA 2 



j^(co-iso) 
^^(co-iso) 

C1+C2 — »(si+S2) 

2 





_ 1A|( C - * s o) 



c 1 +c 2 —i(s 1 +s 2 ) 
2 



n+l 



where we use a notation analogous to (2.18) and write 

(2.26) exp(^iAt) = c\ — is%, exp(^ 2 Ai) = c 2 - «S2, c := c x - c 2 , s := si - s 2 . 

Clearly, the algorithm given above is first order in time. We can get a second order scheme by the Strang 
splitting method, which means that we use Step 1 with time-step At/2, then Step 2 with time-step At, and 
finally integrate Step 1 again with At/2. Our algorithm given above is an 'explicit' and unconditional stable 
scheme. The main costs are DFT and IDFT. 

Lemma 2.1. Our numerical scheme conserves the particle density in the discrete I 2 norm (discrete total 
charge) and the Lorentz gauge. 



Proof: From (|2.17|) and (|2.25|), it is easy to check that the discrete total charge is conserved. From the 

, iAt 



initial conditions and (2.19), we have 





6d t V° 



■it -A =0, ep u 



\ev°-isz-d t A u , p 1 



p 
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From (|2.19|) and (|2.20|), we obtain 



At 2 \£\ 2 
4<5 2 



Tt + l 



?£• A 



n+l 



= I 1 - 
At 

+ T 

sAt 

+ -w 



At 2 \£\ 2 
45 2 



Sd t V n + i£ ■ A 



~\£\ 2 V n + i5£ ■ d t A r 
P n + P r 



26 ? 



Then it is clear that for all n, we have 



Sd t V n + A = 0, ep n = |f| 2 V n - i6£ ■ d t A , p n+1 



iAt 
~25 



- n * n+l 

J +j 



□ 



In order to test the numerical scheme we consider the example of an exact solution for the full MD system, cf. 
[T2"] . In all of the following examples, we take the computational domain fl to be the unit cubic [—0.5, 0.5] 3 . 

Example 2.1 (Exact solution for the MD system). Let us consider the MD system for e = S = 1 with 
initial data 

exp(i£ • x) 



V> (0) (x) = 



(2.27) 



V {0} (x) = V {1) (x) = 0, A (0) (x) = A (1) (x) = 
and external fields given by 

t 2 £ 

a ex S 



x, x = (6,£i+*'6,vTT^-i,o), 



(2.28) 



ye 



2 • 



£ = (2?r, 4tt, 6tt) e 



In this case, there is an exact plane wave solutions for the MD system in the following form, cf. 



(2.29) 



y/2Q. + \Z\*-y/T+W) 

t 2 £ 

-, A(t,x) = — , 



V(t, x ) 



t 



In Figure [T] we see that our method gives a very good agreement with the exact result. 

To test the accuracy of our time-splitting method for the MD system, we did the spatial and temporal 
discretization error tests (see Table [T] and [2]). Table[T]shows the spectral convergence for spatial discretization. 
Table [2] shows the convergence rate for temporal discretization is about 2.0. Here "0 A:r ' At (£, •) is the numerical 
solution for mesh size Ax and time step At, and ip(t, •) is the exact solution given by (2.291. In the following 
also show the charge conservation test (see Table [3| : 

Table 1. Spatial discretization error test: at time t=0.25 under At = 1/1024 (e = S = 1). 



mesh size 


Ax =1/4 Ax = 1/8 Ax = 1/16 Ax = 1/32 




^ A *' A *(f,-)-V(t,0 [2 


8.40E-2 2.68E-3 6.95E-5 5.00E-8 


\W,-)\\p 


convergence order 


4.9 5.3 10.4 
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FIGURE 1. The numerical solutions of example 2.1 Here At — j|g, Ax = The top two 
graphs are the real and imaginary parts of i(^i(t, 0> x 2> 0)|t=i. The bottom two graphs are 
electromagnetic potentials and the norms of ip^, k — 1,2, 3, 4. The solid lines are the exact 
solutions, 'ooo', '***', 'ooo' and are numerical solutions. 



Table 2. Temporal discretization error test: at time t=0.25 under Ax = 1/32 (e = S = 1). 



time step 


At 16 At — 32 At — 64 At — 12g 




^^(V)-V(V) 


I 2 


2.59E-4 5.14E-5 1.29E-5 3.21E-6 


\m,-)\\ P 


convergence order 


2.3 2.0 2.0 



TABLE 3. Charge conservation test: under Ax = 1/32, At = 1/1024 (e = 5 = 1). 



time 


t=0 


t=0.5 t=1.0 




Hp 1.00000000 


0.99999998 0.99999997 



3. The semi-classical regime 



We shall consider in the following the semi-classically scaled MD system (1.14 1. First we shall discuss 
(formal) asymptotic description as e —¥ and then consider some particular numerical test cases. 
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3.1. Formal asymptotic description. To describe the limiting behavior of ip E as e — > we introduce the 
following notations: 

Analogously to the free Dirac operator, the matrix- valued symbol 2? A (£) can be (orthogonally) decomposed 
into 



hid) :=±A A (0 + nx), 



(3.1) 
where 
(3.2) 
with 

(3.3) A A (0 := ^l + ie-A ea; (x)| 2 + F e *(x). 

The corresponding (orthogonal) projectors n A (£) are then given by 



(3.4) 



n A (0 



i 4 ± 



MO 



(2?a(0 - V ea (x)I 4 ) 



Clearly, we obtain the corresponding decomposition of the free Dirac operator (2.5|, (2.7), by setting 
A ex (x) = and V ex (x.) — in the above formulas. Note that /i A (£) is nothing but the classical rela- 
tivistic Hamiltonian (corresponding to positive resp. negative energies) for a particle with momentum £. 
These particles can be interpreted as positrons and electrons, resp., at least in the limit e — > 0, as we shall 
see below. Finally, we also define the relativistic group-velocity by 

(3.5) := V 

The group velocity for free relativistic particles is then (£) = £/Aq(£). 



To do so we assume (well prepared) highly oscillatory 



The semi-classical limit for solution of the weakly nonlinear MD system (1.14) can now be described by means 
of WKB-techniques as given in (see also [! 
initial data for ip e , i.e. 

(3.6) 



V> (0) (x) - ufixyti (x)/£ +u7(x)e' i ^ x )/ £ + 0(e). 



We then expect that ip s (t,yi) can be described in leading order (as e —> 0) by a WKB-approximation of the 
following form 

(3.7) 



V> e (£,x) ~ w+(i,x)e^ (*< x >/ £ +u-(i,x)e^ ( *< x)/e + 0(e). 



Here, the phase functions <\> (£,x) € K, resp. satisfy the electronic or positronic eiconal equation 



(3.8) 



a^frx) + ^(V^^x)) = 0, 0±(O,x) = 0±(x). 



As usual in WKB-analysis we can expect an approximation of the form (3.7) to be valid only locally in time, 
i.e. for \t\ < t c , where t c denotes the time at which the first caustic appears in the solution of (3.8). 



Remark 3.1. We want to stress that the self-consistent fields A e , V e do not enter in (3.8), i.e. the eiconal 
equation is found to be the same as in the linear case. This is due to the weakly nonlinear scaling described 
in the introduction. In particular, i.e. for the Dirac equation without Maxwell coupling, this setting allows 
us to compute the rays of geometrical optics, i.e. the characteristics for (3.8), independently of A e , V £ . 
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It is shown in [23], for the simplified case where A ex (x) 
u (t, x) € C 4 solve a nonlinear first order system, given by 



V ex (x) = 0, that the principal- amplitudes 



1 



(3.9) 

with initial condition 
(3.10) 



(dt + («4(V0 + ) • V)) u+ + - div(w+(V0+))^ 
(d t + (u>q(V<T) ■ v )) u ~ + \ div(uja{V<t>-))u 



^(O.x) :=II±(V^)«j(x) 



t7V + [it] it 4 



The nonlinearity on the r.h.s. of (3.9) is given by 
(3.11) 



V. 



L«j := A- u>± 

where the fields V, .4. are computed self-consistently through 

(3.12) (da - A)V = p°, (d tt - A)A = J°. 
with source terms 

(3.13) p° :=\u+\ 2 + \u-\ 2 , J°:= W +(V0+)| M +| 2 +c o -(V0-)| W -| 2 . 
The polarization of u is henceforth preserved, i.e. 

(3.14) v^fax) = ^(V^)^^^), for all \t\ < t c , 

and we call u + the (semi-classical) electronic amplitude and u~ the (semi-classical) positronic amplitude. 
Note that in this case, i.e. without external fields, we have the simplified relation 

(3.15) 0+(t,x) = -0-(f,x), 



if this holds initially, which we will henceforth assume. The fact that (3.9) conserves the polarization of u , 
is crucial. It allows us to justify the interpretation in terms of electrons and positrons. In other words, the 
WKB-analysis given above shows that the energy-subspaccs, defined via (3.4), remain almost invariant in 
time, i.e. up to error terms of order 0(e). This, so called, adiabatic decoupling phenomena is already known 
from the linear semi-classical scaled Dirac equation [HI [251 US] • However we want to stress the fact that in 
our non- linear setting rigorous proofs so far are only valid locally in time |23j . More precisely, it holds 



(3.16) 



sup 

0<|t|<t c 



0(e), for every < r < t c 



L 2 



On the other hand we want to remark that in the case of the linear Dirac equation, global-in-time results 
are available which also confirm the adiabatic decoupling for all t € K, cf. [231 126j . 



Note that the nonlinearity in (3.9) is purely imaginary. Hence for the densities p 



,±|2 



we find 



(3.17) 



ft^+div (winwy) = o, 



which clearly implies the important property of charge-conservation: 



(3.18) 



/ (p + (t, x) + p (t, x)) dx = const. 



In the case of non- vanishing external fields, i.e. A ex (x) ^ 0, V ex (x.) ^ the system (3.9) becomes much 
more complicated. First w has to be replaced by w A in the above given formulas and second, an additional 
matrix- valued potential has to be added, the, so called, spin-transport term, cf. [91 1251 [2l>]. which mixes 
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the components of each 4- vector (cf. [25] for a broad discussion on this). We shall not go into further 
details here since in our (semi-classical) numerical examples below we shall always assume A ex (x) = and 
V ex (x.) = 0, since we are mainly interested in studying the influence of the self-consistent fields. The only 



exception is Example 3.3 below, where we treat the harmonic oscillator case with V ex (x) = |x| 



Remark 3.2. Strictly speaking, the results obtained in [53] do not include the most general case of non- 
vanishing external fields and mixed initial data, i.e. u^(0, x) =^ 0. Rather, the given results only hold in 
one of the following two (simplified) cases: Either A ex (x.) = V ex (x.) = and m ± (0,x) ^ 0, or: A ex (x.) ^ 0, 
V ex (x.) 0, but then one needs to assume it + (0,x) = 0, or u~(0,x) = 0, respectively. The reason for this 
is that the analysis given in [23] heavily relies on a one-phase WKB-ansatz, which is needed (already on a 
formal level) to control the additional oscillations induced for example through the, so called, Zitterbewegung 
[32] of J e , cf. [H], for more details. 



3.2. Numerical methods for the WKB-system. In order to solve the Hamilton- Jacobi equation (3 



numerically we shall rely on a relaxation method as presented in [17] ■ Then we can solve the system of 



transport equations (3.9| by a time-splitting spectral scheme, similar to the one proposed for the full MD 



system (cf. Section 2.2 1. Using similar notations, suppose that we know the values u ±,n , V n and A r ' 
Step 1. First, we solve the following problem: 

f a t « ± + v-v^) =v{u ± ), 

(d tt -A)V=p°, 
(d tt -A)A=J a , 

by a pseudo-spectral method, where we use the shorthandcd notations 

(3.20) v^) ^^(V^) (gitrS rjiu*) := ^ div^^V^iA 



(3.19) 



First, we take a DFT of (3.191, i.e. 



(3.21) 



dtu^ 1 + i£ ■ v(u ) —fi(u ± ), 
(d tt + |£| 2 )V=p°, 

(d tt + \t;\ 2 )A=J . 

Let us denote by u ±,n , the value of u ± at time t n in Step 1. Then we can find the values of u ±,n+1 , V" +1 , 
and A n+1 by the Crank-Nicolson scheme. After an IDFT, we obtain the values of u ±,n+1 , V" +1 , and A n+1 ■ 

Step 2. It remains to solve the ordinary differential equation 
(3.22) 



with N given by (3.11). Because 7V ± [u] does not change in step 2, we have 

u ± ' n+1 = exp (i/V^MAt) u ± ' n . 
Remark 3.3. We can also use the Strang- splitting method to obtain a second order scheme in time. Again, 



it is easy to see that this algorithm conserves (3.18) 



The solution of the Hamilton- Jacobi equation (3.8) may develop singularities at caustic manifolds, also 
the group velocities ^{S/cj^) and the principal amplitudes become singular. This makes the numerical 



approximation of the transport equations (3.9) a difficult task. Actually, we are not aware of a previous 
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numerical study on such transport equations with caustic type singularities. Our computational experience 



indicates that it is important to conserve the density in the transport problem (3.9), which relies on an 
accurate (high-order) numerical approximation of the terms w^ : (V0 ± ) and div(w (V^*)). However, the 
Hamilton- Jacobi equation is typically solved by a shock capturing type method, which reduces to first order 
at singularities. In order to get a better numerical approximation, we still use a shock capturing method, 



namely the relaxation scheme developed in |17) . spatially for the Hamilton-Jacobi equation (3.8|, but use the 



fourth order Runge-Kutta method temporally. For the transport problem (3.9) we found that the pseudo- 
spectral method behaves better than finite difference schemes. 

3.3. Numerical examples in the semi-classical regime. In all of the following examples we shall assume 
for simplicity 

(3.23) V (0) (x) = V {1) (x) = 0, A (0) (x) = A (1) (x) = 0, 

since different, i.e. non-zero, initial conditions would only add to the homogeneous solution of the corre- 
sponding wave equation. 

Remark 3.4. Remark that in the following numerical examples <fii has to be chosen such that it satisfies 
the periodic boundary conditions. 



Example 3.1 (Self-consistent steady state). Consider the system (1.14) with initial condition 
(3.24) 



t=o 



X exp 



ixr 

4d 2 



X = (1,0,0,0), d=l/16, 



and zero external potentials, i.e. A£ x (x) = V ex (x.) = 0. This example models a wave packet with initial 
width d and zero initial speed, propagating only under its self-interaction. Note that in this case 0j^(x) = 



and w + (0,x) is simply given by (3.24), whereas it - (0,x) = 0, hence w~(i,x) = 0, for t > 0. First, we choose 
e = 10~ 2 and compare the solution of the full MD system with the numerical solution obtained by solving 
the asymptotic WKB-system (3.8), (3.9). From Figure |2]we see that the two numerical solutions agree very 
well for such a small e. In particular the creation of positrons in the full MD system is small, i.e. O(e) as 
one expects from the semi-classical analysis. This is clearly visible in cf. Figure [3j which shows that the 
projectors (V0) are indeed good approximations of nj(— ieV) for e is small. However for e = 1 this is no 
longer true. Furthermore, because in this case the WKB-phase is found to be simply given by 4> + (t, x) = —t, 
we thus have V(/> + = and V • ujq = 0, and hence the transport equation (3.9) simplifies to 



dtu + + iVu^ 



0. 



which implies |u + (i, x)| 2 to be constant. In this particular case, we can use a very coarse mesh to get 
satisfactory results (cf. Table [4]). Remark that the results in Table [4] also illustrate the validity of (3.16). 
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|2 -- j ^ -,±r+ v ^ c # ± /£| 



IV^x)^^ and E± u (*> x ) 



at i = 0.25 




V e {t,x)\ X3 = and VCt.x)!^ at i = 0.25 




k e (t,x)| X3=0 and J]± w (^ x )e 



at t = 0.5 




V E (t,x)\ X3=0 and F(t,x)| X3=0 at t = 0.5 



FIGURE 2. Numerical results for example |3.1| The left column shows the graphs of the solu- 
tion of the MD system, the right column shows the graphs of the solution of the asymptotic 
problem. Here e = 0.01, At = rko, A.x = 1 



128 ' 



32 • 
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e = 1.0, |n+(-i£V)V E (t,x)| X3= o| 2 and |n+(V0)V £ (t, x)U 3=0 | 2 



FIGURE 3. Numerical results of the densities of electron/positron projectors for example 3. 1 
The left column is |IIq (— zeV)i/' 6 (i, x) | a . 3=0 1 2 , the right column is |IIq (V<^)V> e (t,x)| X3= o| 
Here t = 0.25, Ax = 1/32, At = 1/128. 
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Table 4. Difference between the asymptotic solution and the full MD system for example 
31] (At = 1/128, Ax = 1/32): 



sup 

0<t<0.25 



0.0001 0.001 0.01 
3.20E-3 3.34E-2 2.98E-1 



L 2 (0)«iC 4 



sup 

0<t<0.25 



4.90E-3 5.01E-2 4.40E-1 



Example 3.2 (Purely self-consistent motion). In this example, again zero external fields are assumed, 
but we modify the initial condition for as follows: 

.M_ x )' 



(3.25) 



n= 



X(x) exp 



1/16, 



where the phase function describing the e-oscillations is given by 

1 



(3.26) 



<Mx) 



40 



(1 + cos 27rai)(l 4- cos 2-KX2) 



and we choose the initial amplitude such that iI(j~(V</>j(x))x(x) = x(x), i.e. 
' ff(x)+g(x) 

2(^TT^2_i 



(3.27) 



X(x) 



&(x)(fr(x)-H&(x)) Q gi(x)+i&(x) 

2(7iT|f-i) ' ' 2 



e = V0,(x). 



As in the above example we thus have u~(t,x) = 0. Note that for <f>i = 0, (3.25) reduces to (3.24). The 
numerical solution of the eiconal equation (3.8) [17 indicates a kink-type singularity in the phase of our 
asymptotic description at about t ~ 0.56, cf. Figure [4] Hence the asymptotic WKB-type approximation for 
the spinor field is no longer correct for t > 0.56, 

The numerical results for both the MD system and the semi-classical limit for e = 0.01 are given in Figure 
[5] Table [5] attempts to show the validity of (3.16). Compared to Table |4j the difference between two 
systems is somewhat larger than 0(e). Our experience indicates that this has to do with the numerical 
difficulties mentioned before and with the fact that discretization errors "pollute" the solution of the semi- 
classical system as time evolves, preventing a more accurate comparison at later time. Due to our computing 
capacity, we are unable to conduct more refined calculation, which would have provided a better justification 
of the ansatz (3.16) for this problem. For the same problem, we also present the numerical solutions of the 
Maxwell-Dirac system at later time in Figures [6j We also present a numerical simulation of the case e = 1.0, 
i.e. away from the semi-classical regime, see Figure [8] From the plots it becomes clear that the "exact" 
spinor field and the solution of the asymptotic WKB-problem are qualitatively "close" for small values of e 
and before caustics, while they are even qualitatively different away from the semi-classical regime. 
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Figure 4. The graph of the phase 4>(t,x)\ Xa= o at t = 0.5625 for example 3.2 It shows the 
phase becomes singular at the tip. 



Table 5. Difference between the asymptotic solution and the full MD system for example 
(At = 1/128, Ax = 1/64): 



0.01 0.1 



sup 

0<t<0.125 



sup 

0<t<0.125 



<4 h 



L 2 (n)igic 4 



L°°(fi)<»C 4 



0.196 0.926 



0.115 0.646 



Example 3.3 (Harmonic oscillator). Finally, we take A ex (x) = and include a confining electric poten- 
tial of harmonic oscillator type, i.e. V ex (x) = |x| 2 . Hence ± satisfies 



(3.28) 



flt^fox) + v/|V0±| 2 + l + |x| 2 = 0, 0±(O,x) = ^(x), 



which implies = w o (£) m this case. Due to the presence of the external potential, the semi-classical 



transport equations (3.9) have to be generalized by including a spin-transport term, cf. [26] . which however 



only enters in the phase of u . Thus the conservation law for the densities p ± is the same as in (3.171 



Let us consider the system (1.141 with initial condition 

(x 1 -0.1) 2 + (x 2 +0.1) 2 + a;2 



(3.29) 



^ \ t=0 = X exp 



4d 2 



X = (1,0,0,0), d= 1/16, 



In this case we choose e = 10 2 , At = 1/32, Ax = 1/32. The numerical results are give in Figure[9j We see 
that the wave packet moves in circles due to its interaction with the harmonic potential. 

Remark 3.5. In analogy to the spectral-splitting method for the Schrodinger equation analyzed in pQ, we 
find that Aa; = 0(e) and At = 0(1), as e —> 0, is sufficient to guarantee well-approximated observable of 
the MD system. A more refined grid in temporal direction is necessary to obtain a good approximation for 
the reps, components of the spinor field itself, typically At — 0(e 2 ) is needed. 
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\ip e (t,x.)\ X3 = \ and J2± u (*> x ) 



x 3 =0 



at t = 0.25 




-0.5 -05 -0.5 -0.5 

V £ (t,x)\ X3=0 and V(t,x)\ X3=0 at t = 0.25 




|V/(t,x)U 3= o| 2 and £. u^xje**/ 6 



at i = 0.375 




V e (t,x)| a ,=o and V(t,x)U,= at t = 0.375 



Figure 5. Numerical results for example 3.2 The left column shows the graphs of the solu- 
tion of the MD system, the right column shows the graphs of the solution of the asymptotic 
problem. Here e = 0.01, At = y|g, Ax = 
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I U 3=0 1 and V e (t,x)| X3=0 at t = 0.53 




Re(V>i(£,x))| X3= o and Im(^f (i,x))| X3=Q at i = 0.53. 




|^ £ (i,x)| X3=0 | and ^(t.xJI^o at t = 0.625 




Re(^i(*,x))U 3 =o and Im(^f (i,x))| X3=0 at i = 0.625. 



Figure 6. Numerical results of the MD system for example 3.2 Here e — 0.01, At = ^ 



128 ' 



Ax 32 . 
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i = 0.25, |noHeV)V e (t,x)| X 3 = o| 2 and |ll^ (V0)^(i, x)L =0 | 2 




* = 0.25, |n+H e V)V> e (t,x)| X3= o| 2 and |n+(V^ e (t,x)| X3=0 | 2 




t = 0.75, I (-ieV) V E (*, x) U 3=0 1 and |n^(V^ £ (t,x)| X3 = | 





1 

08 






06 






04 

02 





f = 0.75, |n+HeV)V e (t,x)| X3= o| 2 and |n+(V0)^(t,x)| X3=o | 2 



Figure 7. Numerical results of the densities of electron/positron projectors for example 



3.2 



The left column is |II (— feV)^> e (t, x)| X3= o , the right column is |±I (V0)?/> E (t)| X3= o 
The graphs show that the matrices n^"(V</>) do not mimic n^"(— ieV) after the caustic point. 
Here e = 0.01, Ax = 1/32, At = 1/128. 
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-0.5 -0.5 -0,5 -0.5 -0.5 -0.5 

Re(<^f(f,x)) 1^=0, Im(V>i(t,x))| X3=0 and V(t,x)\ X3=0 at t = 0.25. 

\fy - if f 



Re(^f(t,x)) 1^=0, lm(-0f(i,x)) | K 3 = o and V(t,x)\ X3=0 at i = 0.5. 



Re(^i(*,x)) 1^3=0, Jm(il>l(t,x))\ Xa= Q and V(t,x)\ X3=0 at t = 0.625. 



Figure 8. Numerical results of the MD system for example 3.2 Here s = 1.0, At = jj^, 
Ax = 



32 • 



4. The non-relativistic regime 



Finally we shall also consider the non-relativistic regime for (1.18) as S — > 0. Again we shall first describe 
the formal asymptotics and then discuss numerical examples. 

4.1. Formal description of the asymptotic problem. To describe the non-relativistic limit of the MD 
system we first define two pseudo-differential operators 11*, (Z?) via their symbols 



(4.1) 



e/pV 



where Ao(£), are given by (2.6), (2.4). We then define the (non-relativistic) electronic and the (non- 



relativistic) positronic component ip^, ipp by 

(4.2) ^(t,x) :=e it /« a n2(D)^(i,x), ^(i,x) := e^^D)^, x), 

where V'' 5 solves the non-relativistically scaled MD system (1.18). Note the difference in sign of the phase- 



factors. This corresponds to subtracting the rest energy, which is positive for electrons but negative for 
positrons, cf. (HJ [H 120) . The above given definition of electronic/positronic wave functions should not 
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-U.5 -0.5 -0.5 -0.5 -0.5 -0.5 

t=3 t=4 t=5 

1 .5 1 .5 

^^^^^°' 5 ^^^^^°' 5 ^^^^^ a 

-0.5 -0.5 -0.5 -0.5 -0.5 -0.5 

t=6 t=7 t=8 

5 1.5 1.5 



Figure 9. Numerical results of the density at different time for example |3.3| Here e = 0.01, 
Ax = 1/32, At = 1/32. 



be confused with the one obtained in the semi-classical regime, since both definitions are adapted to the 
particular scaling of the resp. system under consideration. We remark that up to now there is no satisfactory 



interpretation in terms of electrons and positrons for the solution of the full MD system (1.1), (1.2). Indeed 



there is no such interpretation even for the linear Dirac equation with external fields, see e.g. 
Remark 4.1. It is easy to see that the formal limit 8 — > of the operators 11*/ (D) yields, 



(4.3) 



n 



% o^ 
. o 0, 



IT 



'0 
,0 I, 



This explains the interpretation of electrons (resp. positrons) as the upper (resp. lower) components of the 
4- vector ip s for small values of 5, cf. [2"2"] . 

It is then shown in [5] (see also [B] for easier accessible proofs in the linear case) that 



(4.4) 



5->0 



^ /p (t,x)^> e/p (<,x), inff 1 ^ 3 ) 
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(4.5) 



< 



idt<p p 



v ex ) Vp , 



where ipeifp solve the mixed electronic/positronic Schrddinger-Poisson system: 

A 
2~ 
A 

j<Pp + (V 
-AV = \<p p \ 2 + \p e \ 2 , 

In contrast to the asymptotic problem obtained in the semi-classical limit, this system is globally well posed. 
The appearance of the Poisson equation can be motivated by performing a naive Hilbert expansion in the 
self-consistent fields, cf. [H], i.e. 



(4.6) 



V 



V + SV + 0(S 2 ), A = A + SA + 0(8 2 ). 



Plugging this into (1.18), comparing equal powers in 5, and having in mind that J 

as S 



0(1) [6 gives (4.5) 



In [5] the electric potential is proved to converge in iJ^ll 3 ) as 6 —> 0, whereas the convergence of the 
magnetic fields is not studied in detail. Indeed, it is shown in [8] that if one only aims for a derivation of the 
Schrodinger-Poisson system, one can even allow for initial data A 5 (0,x), d t A s (0,x) which do not converge 
as 5 -> 0. 

Remark 4.2. If we would, in addition, consider terms of order 0(6) too, we (formally) would obtain a 
Pauli equation for <p e / p , including the matrix-valued magnetic field term ^2<7kBk, i-e. the, so called, Pauli- 
Poiswell system, cf. [6l [19]. Moreover we remark that the authors in [8 considered the MD system in 



Coulomb gauge, i.e. divA = 0, instead of the Lorentz gauge condition imposed in this work (1.7). The 



reason is rather technical and it is not clear yet if a generalization of their work to the Lorentz gauged system 
is possible. 

As before we shall use a time-splitting spectral method [2] to solve the coupled system of Schrodinger-Poisson 



equations (4.5) 



Step 1. First, we solve the following problem: 

id t <p e 



9<- 



\9e\\ 



(4.8) 



\-AV = \<p p \ 2 
Step 2. Then we solve the coupled equations 

id t 9 e = (v + v ex )v e , 

id t <p p = (V + V ex )<p p , 

In step 1, we again use the pseudo-spectral method. In step 2, we can get the exact solution for this linear 
ODE system in time, since |</5p| 2 and |<p e | 2 j resp., are kept invariant by step 2. 

Remark 4.3. Let us fix e — 1 and consider S in the algorithm given in section |2.2| Based on the 



expansion of (2.17)-(2.20), we obtain 

=exp(A(t -*„))*" + 0(<J), 



(4.9) 

where in the limit S — > the matrix A € C 4x4 simplifies to 

A = diag[A,A, -A, —A], A = -i(S~ 2 



ICIV2). 
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We also have 
(4.10) 
and 
(4.11) 



|£| 2 (y n + v n+1 



|$n|2 + |$n+l|2 + 0(*) 



. n+1 



0(S), 



because a fe $ n+1 ) = 0(6). If we denote the upper (resp. lower) components of the 4-vector ^ by \& e 

(resp. ^ p ), we obtain 

M\ 2 



(4.12) 



ft (e ±lt 



/p 



/p 



0(<5), 



and from (2.25), we find 
(4.13) 



n+l 



exp (-iVAt) $ 



n+l 



0{8). 



Combining the equations ( 4.10 1— ( 4.13 ) , we conclude that the numerical solutions of our algorithm, given in 



section [2T2| uniformly converge to the numerical solutions of the above algorithm. This analysis, previously 
done for a time-splitting spectral method for the Zakharov system [IS] , shows that one can choose h, At 
independent of S. 

4.2. Numerical examples for the non-relativistic regime. 



Example 4.1 (Purely self-consistent motion II). Here we consider the MD system (1.18) in a unit 
cubic with periodic boundary conditions, zero external fields, and initial data 

2 - 



(4.14) 



^(x)| t=0 = V (0) (x) = X exp 



1*1 
Ad 2 



x = (1,1,1,1), d = 



l 

16' 



A^°>H^°>| 2 , ^ (1) (x) = 0, 

A4 0) = (^°>,a fe V (0) )c4, A«(x) = 0. 



Note that the above choice of initial data for V and Ak is done to avoid initial layers. The impact of this 
choice on the numerical resolution, i.e. the mesh strategy etc., is analogous to the Zakharov system discussed 
in [TS]. We also consider the Schodinger-Poisson problem (|4.5|) with the initial data 

(4.15) 



<^(t,x)| t=0 = nf(£>)^°>(x), ¥> p (t,x)| t=0 = Tl s p (D)^°\x), 



We compare the solution of the MD system with the (coupled) Schodinger-Poisson problem, cf. Figure 10 



and Figure [12] Table [6] Fi gure and Figure [12] illustrates the validity of (| 4.4[ ) . The Figures |10[]13| also 
show that |A 5 | = 0{S)\V s \, as 5 -> 0. 



Table 6. Convergence test for example 4.1 (here At = 1/128, Ax = 1/64) 



5 


0.01 


0.1 


1.0 


SUp - f + \l/lp - fp\ 2 
0<t<l/4 


0.101 


0.345 


2.407 
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|^(i,x)| s 



1^3=0 



Re^^^x))! Im^ft.x))! and V*(i,x)| for J = 1.0. 




|V/(t,x)| 2 | Re^^x))! Im «i(t,x))| and ^(t,x)| for <S - 0.01. 




(|<^ e (t,x)| 2 + |^(i,x)| 2 ) | x3=0 , Re(v3 e ,i(t,x))| a;3=0 , Im(^ e) i(t,x))| !e8=0 and V(t,x)| xa 



Figure 10. Numerical results for example |4.1| at t=0.5. The first row is the solution of the 
MD system with 5 = 1.0, whereas the second row is the solution of the MD system with 
S = 0.01, the third line is the solution of the Schrodinger-Poisson system. 




A?(t,x)| A|(t,x)| x3=0 , A«(t,x)| for 5 = 1.0. 




Af(t,x)| 4(i,x)l =Q , 4(t,x)l Q for 5 = 0.01. 



Figure 11. Numerical results of the magnetic fields for example 4.1 at t=0.5. The first 
row is the solution of the MD system with 5 — 1.0, whereas the second row is the solution 
of the MD system with 5 = 0.01. 
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|^(i,x)| 2 | x3=0 , Re^li^x))!^, Im^t.x))!^ and Vfox)]^ for 5 = 1.0. 




|V/(t,x)| 2 | Re^^x))! Im^^x))! and V s (t,x)\ for 5 - 0.01. 




(\ip e (t,x.)\ 2 + |<^(i,x)| 2 ) | a3=0 , Re(v3 e ,i(t,x))| :E3=0 , Im ((pe,i(t, x )) L 3=0 and ^> x )L=cr 



Figure 12. Numerical results for example |4.1| at t=1.0. The first row is the solution of the 
MD system with S = 1.0, the second row is the solution of the MD system with S = 0.01, 
and the third row is the solution of the Schrodinger-Poisson problem. 




Af(t,x)\ Q , Af(t,x)| 4(*,*)L =0 ^ S = 1.0. 




A?(t,x)| ^(* s x)| Bs=0J ^(i,x)| for J = 0.01. 



Figure 13. Numerical results of the magnetic fields for example 4.1 at t=1.0. The first 
row is the solution of the MD system with S = 1.0 and the second row is the solution of the 
MD system with S = 0.01. 
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|V 5 (i,x)| 2 | x3=0 , (|<^(i,x)| 2 + Mt,x)| 2 ) I V*(t,x)| and V(t,x)\ at t=0.5. 




|V> 5 (i,x)| 2 | (|^(t,x)| 2 + |^(t,x)| 2 ) I V*(*,x)| and V(t,x)\ at t=1.0 




l^(*,x)| 2 | (|^ e (t,x)| 2 + |<ft,(t,x)| 2 ) | V s (t,x)\ and V(t,x)| at t=1.5 




l^(^ x )| 2 L 3=0 ' (IVe(t,x)| 2 + |<^(t,x)| 2 ) | K3=Q , y 5 (t,x)| x3=0 and V(t,x)| X3=0 at t=2.0. 



Figure 14. Numerical results of the density for example |4.2| The first and third column are 
\ip s (t, x)| 2 | and V^(i,x)| for MD system, respectively. The second and fourth col- 
umn are (|y e (t, x)| 2 + |<p p (i,x)| 2 ) | _ Q and V(t,x)| _ Q for Schrodinger-Poisson equation, 
respectively. Here 6 = 0.01, Ax = 1/64, At = 1/128. 



Example 4.2 (Harmonic oscillator II). Finally, we choose A ea: (x) = but include a confining electric 
potential of harmonic oscillator type, i.e. V ex (x) = C|x| 2 . To compete with the effect of the diffusion term 



Aip s , we choose the large constant C = 100. Let us consider the system (1.18) with initial condition 
(4.16) 



ft* I I (zi-0.1) 2 + (x 2 + 0.1) 2 + z 2 \ , . ... 

^ \t=o = X ex P ( ) ' X = (1, 0, 1, 0), d = 1/16, 



In this case we choose 8 = 10 , At = 1/128, Ax — 1/64. The numerical results are shown in Figure 14 We 



see that the wave packet moves in circles due to its interaction with the harmonic potential and the diffusion 
term A-0 5 . Note that agreement with the non-relativistic results is very good also for this test. 
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5. Conclusion 

In this work, we presented a time-splitting spectral scheme for the MD system and similar time-splitting 
methods for the corresponding asymptotic problems in the (weakly nonlinear) semi-classical and in the non- 
rclativistic regime. The proposed scheme conserves the Lorentz gauge condition, is unconditionally stable 
and highly efficient as our numerical examples show. In particular, we presented numerical studies for the 
creation of positronic modes in the semi-classical regime as well as numerical evidence for the smallness 
of the magnetic fields in the considered non-relativistic scaling. A distinct feature of our time-splitting 
spectral method, not shared by previous methods (using the time-splitting spectral approach) , is that in the 
non-relativistic limit, the scheme exhibits a uniform convergence in the small parameter 5. 

We finally remark that there are several open questions that deserve further exploration. For example, it 
would be an interesting project to derive a better numerical method for the system of eiconal and transport 
equations, describing the semi-classical limit, which consequently would allow for a more accurate comparison 
between the solution of the MD system and limiting WKB-description. A second step then should be 
the numerical study of the semi-classical MD equations with stronger nonlinearities, in particular 0(1)- 
nonlinearities, a so far completely open problem, even from an analytical point of view. 
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